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Thermal fluctuations of microtubules (MTs) and other cytoskeletal filaments govern to a great 
extent the complex rheological properties of the cytoskeleton in eukaryotic cells. In recent years, 
much effort has been put into capturing the dynamics of these fluctuations by means of analytical and 
numerical models. These attempts have been very successful for, but also remain limited to, isotropic 
polymers. To correctly interpret experimental work on (strongly) anisotropic semiflexible polymers, 
there is a need for a numerical modelling tool that accurately captures the dynamics of polymers with 
anisotropic material properties. In the current study, we present a finite element (FE) framework for 
simulating the thermal dynamics of a single anisotropic semiflexible polymer. First, we demonstrate 
the accuracy of our framework by comparison of the simulated mean square displacement (MSD) 
of the end-to-end distance with analytical predictions based on the worm-like chain model. Then, 
we implement a transversely isotropic material model, characteristic for biopolymers such as MTs, 
and study the persistence length for various ratios between the longitudinal shear modulus, Gi2^ 
and corresponding Young's modulus, Ei. Finally, we put our findings in context by addressing 
a recent experimental work on grafted transversely isotropic MTs. In that research, a simplified 
static mechanical model was used to deduce a very high level of MT anisotropy to explain the 
observation that the persistence length of grafted MTs increases as contour length increases. We 
show, by means of our FE framework, that the anisotropic properties cannot account for the reported 
length-dependent persistence length. 



I. INTRODUCTION 

The characteristic thermal fluctuations of biopolymers 
that are submerged in a viscous solvent have been ex- 
tensively studied for over a decade [IHl], mainly due to 
their important role in many biological processes [3 |6] . 
In eukaryotic cells for example, the thermal undulations 
of biopolymers are essential for the functioning of the dy- 
namic network that these biopolymers constitute [71 18] . 
This network is known as the cytoskeleton and mostly 
consists of the biopolymers F-actin, intermediate fila- 
ments, and microtubules, which are linked together by 
accessory proteins. The cytoskeleton gives the eukary- 
otic cell its mechanical properties, allows the cell to re- 
sist external stresses and plays a leading role in active 
force generation, cell mitosis and intracellular transport 
of vesicles and organelles [H [9] . A good understanding 
of the mechanisms underlying these divergent roles of 
the cytoskeleton, requires a complete description of the 
mechanics and thermal dynamics of all of its structural 
components [TT]. 

The thermal dynamics of a submerged filament is of- 
ten discussed in terms of its persistence length, Ip. If 
the contour length, Lc, of the polymer is much greater 
than its persistence length, the dynamics of the polymer 
is dominated by thermal bending (flexible limit). Con- 
versely, if the contour length is much smaller than the 
persistence length, the polymer dynamics will be domi- 
nated by elastic bending (stiff limit). Many biopolymers 
have a contour length close to their persistence length 
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and are therefore called semiflexible. Formally, the per- 
sistence length is the characteristic length scale of the arc 
length for which two tangent angles along the polymer's 
backbone 6{0) and 6{s) start to become uncorrelated 
In two-dimensions, this translates to [3] : 

(cos(i9(5) - 6>(0))) = exp{-s/2£p) . (1) 

By means of the equipartition theorem, the persistence 
length can be related to the bending rigidity of the poly- 
mer, hz, and the external thermal energy, /c^T, as [5 

£p = K/ksT . (2) 

For homogeneous isotropic filaments. Hi equals the prod- 
uct (EI) of the Young's modulus, E^ and second mo- 
ment of area, /. However, the effective bending stiffness 
of anisotropic polymers is nontrivial and cannot be ex- 
pressed in terms of individual material properties [iT. 

Microtubules (MTs) are among the cytoskeletal 
biopolymers that have attracted particular interest dur- 
ing the recent years, mainly due to their key role in cell di- 
vision [12] . They exhibit an intrinsic anisotropic material 
behaviour [2 [3l [13] . Underlying this anisotropic material 
behaviour is their complex molecular structure: MTs are 
hollow tubes with an outer diameter of 25 nm that con- 
sist of parallel aligned protofilaments of a- and /3-tubulin 
heterodimers (Figure [T]). The longitudinal Young's mod- 
ulus is determined by the strength of the head-tail a/3-af3 
tubulin bonds, whereas the shear modulus is determined 
by the much weaker inter-protofilament bonds p!4HT6] . 
All- atom computer simulations [17 and in vitro mechan- 
ical testing of MTs [13 report anisotropic properties with 
one order difference between the longitudinal Young's 
modulus and corresponding shear modulus. 
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FIG. 1. A schematic image of the microscopic architecture of 
microtubules. Microtubules are straight, hollow filaments and 
consist out of, on average, 13 parallel aligned protofilaments. 



Mainly due to the unfathomable relationship between 
anisotropic material properties of a polymer and its per- 
sistence length, it has been difficult to interpret recent 
experiments on grafted MTs submerged in a viscous fluid. 
In those experiments, a hidden dependency of the persis- 
tence length of grafted MTs on their contour length was 
observed [21 118] . Based on a simple Timoshenko beam 
model, that describes macroscopic elastic beams with 
shear contributions [19 , and by assuming a static tip 
load and small displacements, anomalous MT material 
properties of a six orders lower longitudinal shear mod- 
ulus than corresponding Young's modulus were found. 
Various numerical modelling studies indeed confirmed 
contour length-dependency of the persistence length, but 
these studies are likewise limited to static loads and small 
displacements [20'-'22l. In this respect, there is a lack of 
an intuitive computational modelling tool that allows for 
the study of the thermal dynamics of anisotropic poly- 
mers under distributed thermal force application and fi- 
nite displacements. 

Conventionally, polymers have been modelled numer- 
ically as a series of interconnected beads and rods or 
springs. This approach has its limitations, such as the 
lack of careful mathematical analysis, unavoidable arti- 
ficial constraints and the need for (expensive) explicit 
time integration schemes [23 . Furthermore, bead mod- 
els are limited to very slender and isotropic filaments, 
which makes these models less suitable for modelling 
the dynamics of relatively thick and anisotropic poly- 
mers, such as MTs. Driven by these limitations, a new 
method based on the finite element (FE) method was 
recently proposed |23]. In such an FE model, the mi- 
crostructural complexity of a polymer is approximated 



by a continuum mechanical model, which is discretized 
into finite elements. The time-dependent partial differ- 
ential equations (PDE's), resulting from the interaction 
between the polymer and its fluid environment, are solved 
over each element using implicit time integration. It can 
be proven mathematically that the solution of the FE 
method converges to the solution of the analytical PDE 
for infinitely small elements and time steps [23 . Addi- 
tionally, FE modelling has already been used for many 
decades in other research fields, such as computational 
engineering. This has resulted in versatile user-friendly 
software packages, allowing straightforward implemen- 
tation of advanced material models. Furthermore, the 
FE method has been shown to be up to thousand times 
faster than traditional bead-rod models based on explicit 
time integration schemes [23^. The sound mathemati- 
cal formulation, easy implementation of complex mate- 
rial models and favourable computational costs makes 
the FE method an accurate and intuitive modelling tool 
for capturing the thermal dynamics of single semiflexible 
polymers. 

In the current study, we develop a finite element frame- 
work, based on a commercial solver (Abaqus/ standard, 
SIMULIA), to model semiflexible polymers in thermal 
equilibrium with their viscous fluid environment. We first 
demonstrate the accuracy of this technique by comparing 
the simulated time evolution of the thermal fluctuations 
of freely floating and hinged isotropic filaments with an- 
alytical predictions based on the worm-like chain (WLC) 
model. We then change the boundary conditions to sim- 
ulate a grafted MT and use the mean square of the trans- 
verse displacement of its tip as a measure for its persis- 
tence length. Again, we show good correspondence with 
the theoretical predictions based on the worm-like chain 
(WLC) model. Finally, we implement highly anisotropic 
material properties and study the persistence length of 
MT of various lengths. We relate our findings to recent 
experimental research [2^. 



II. METHODOLOGY 
Introduction 

In this section, we describe the numerical framework 
for capturing the thermal dynamics of single anisotropic 
semiflexible polymers. This numerical framework is 
built around a commercial FE solver (Abaqus/Standard, 
SIMULIA) and is largely based on the framework proposed 
by Cyron and Wall [23 . The finite element approach 
is, as argued by Cyron and Wall in [23l [24], an accu- 
rate, efficient and intuitive way of modelling the ther- 
mal dynamics of semiflexible polymers. In the current 
section, we first explain the workflow of the framework 
centred around its three main stages (pre-processing, 
solving and post-processing). Then, to confirm our ap- 
proach, we pay particular attention to the validation of 
the framework. Based on the time evolution of the mean 
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square displacement (MSD) of the end-to-end distance of 
simulated isotropic MTs, we show good correspondence 
between simulation output and the analytical solution 
based on the WLC-model. We do this for various pa- 
rameters and boundary conditions. Finally, we intro- 
duce an anisotropic material model and adapt the vali- 
dated framework to address a recent experimental study 
on grafted anisotropic MTs [2 . 
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A. The finite element framework 

From here on, we coarse-grain the exact atomic archi- 
tecture of an MT by approximating it as a hollow slender 
continuum structure, as is common practice in Brownian 
dynamics modelling [2^ |5j . Let us formulate a force equi- 
librium per unit length of internal and external forces to 
which the slender continuum is subjected |23], 



finert(u) +fint(u,U,x) = fext(x). 



(3) 



In Eq.([3|, the internal forces, fint, include bending (elas- 
tic) forces, fei(u, x), and damping forces, fvisc(u, u, x). 
Since we do not consider any deterministic external 
forces, the right hand side of Eq.(|3]) only includes stochas- 
tic forces, fstoch(x). Furthermore, from dimensional anal- 
ysis it can be shown that on the scale of a polymeric 
system the inertial forces are negligable compared to the 
internal friction, elastic and stochastic forces [25 . There- 
fore, from here on we will not explicitly write the inertia 
force term. We thus rewrite Eq.([3| as [23] 



fel(u, X) + fvisc(u, U, X) = fstoch(x). 



(4) 



Eq.(|4| is a nonlinear partial differential equation that 
cannot be solved analytically. We therefore resort to EE 
modelling for finding a solution for the displacement vec- 
tor, u. In order to efficiently perform the simulations, 
a modelling framework was built around a commercial 
EE solver (Abaqus/Standard, SIMULIA), based on a 
Fortran user-subroutine, MATLAB and Python-scripting 
(Table [l|. 



TABLE I. Overview of the software that was used in the finite 
element framework 



Software 


Version Manufacturer 


Function 


Python 


2.6.2 Beopen.com 


pre / post-processor 


Abaqus / Standard 


6.10-1 SIMULIA 


solver 


Fortran 


11 Intel 


solver 


MATLAB 


R2011a Mathworks 


post-processor 



1 . Pre-processing 

Let us consider a polymer confined to movement in 
the x-y plane (2-D). At tsim = 0, the backbone of the 



FIG. 2. Initially the filament is aligned with the x-axis. Move- 
ment is confined to the x-y plane and the biopolymer (top) is 
discretized into a one-dimensional model of 20 elements (bot- 
tom). 



polymer is aligned with the x-axis, which can be consid- 
ered the stress-free or reference configuration [26]. The 
contour length of the polymer, Lc, is a variable param- 
eter and is specified by the user for each ensemble. The 
polymer is then discretised into 20 elements of length 
le = Lc/20 (Figure|2|. Cyron and Wall (2009) showed al- 
ready excellent results for discretization into 20 elements 
to capture the thermal dynamics of semifiexible poly- 
mers of similar properties as in the current study, and 
only minor improvements were reported a for more accu- 
rate discretization [23 . For discretization we used one- 
dimensional Timoshenko beam elements. This specific 
type of elements was chosen due to their ability to accu- 
rately describe shearing, even for relatively thick beams 
and low shear modulus [27 and their favourable compu- 
tational efficiency [26]. Additionally, several studies have 
shown that one-dimensional Timoshenko beam elements 
can accurately capture microtubule mechanics and are 
preferable over orthotropic shell models [28l [29]. The 
cross section of the filament is dependent on the type of 
polymer that is simulated and is presumed fixed through- 
out different ensembles [30 . Since we simulated MTs, 
which are of tubular shape, we assumed a hollow circular 
shaped cross-section, of which the dimensions are given in 
Table [llj We used one-dimensional continuum elements, 
thus the MT cross-section is only implicitly reffected by 
the second moment of area, / [26l: 



4 4 
1 = 7T , 



(5) 



where Vq is the outer radius and is the inner radius of 
the polymer. First, for validation of the FE framework, 
we implemented an isotropic material model by setting 
the shear modulus to Gf§ = Ei/2{1 + u). See Table |ll| 
for a listed overview of all the geometric and material 
properties that are used in this study. 

Stochastic forces are applied on each individual node 
of the discretized polymer by means of a user-defined dis- 
tributed load subroutine (DLOAD). At the beginning of 
each simulation time step AtpE the DLOAD subroutine is 
called, upon which it outputs an array of random forces. 
These discrete forces are the equivalent of the noise corre- 
lations, fixed through the fiuctuation-dissipation theorem 
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TABLE II. Values that were used for validation of the mod- 
elling framework 



Parameter 




Value 


Ref 


WU-Uci icLLllU-O 


To 


1 oc; V 1 V\~^ nm 

L.ZiO A ±U JJ-lii 


M 


Inner radius 


Ti 


7.5 X 10"^ pm 


u 


Poisson's ratio 


V 


0.3 


im 


Density 


P 


1.0 X 10^ mgpm"^ 


m 


Young's modulus 


E 


1.3 X lO^Nyim-^ 


m 


Fluid viscocity 


V 


1 X 10"^pNsyim"^ 


m 


Thermal energy 




4.045 X 10"^pNyim 


m 



[23^, ^25^ , that determine the continuous stochastic force 
density fstoch- These random nodal forces are sampled 
from a Gaussian distribution, which is defined by its first 
and second moments, 



/i = 

A^stoch^e ' 



(6) 



respectively, where we introduced the thermal energy of 
the solvent, /c^T, the drag-coefficient, and the variable 
A^stoch, which represents the 'refresh rate' of the stochas- 
tic forces. The FE simulation step size, AtpE, has to be 
smaller than Atgtoch in order to obtain proper FE con- 
vergence. In the current framework, we have consistently 
used a FE step size AtpE ^ (0-lAtstoch)- 

Similar to earlier work [23^ , we estimated the homoge- 
neous isotropic drag coefficient, according to the for- 
mula for a rigid cylinder in a homogeneous flow [25 as. 



C = 47r7?/ln(Vrf) , 



(7) 



where 77 is the fluid viscosity and d is the diameter of the 
filament [26 . The logarithmic term is a correction factor 
that compensates for the neglect of hydrodynamic inter- 
actions between distant segments of the polymer [25] |33] . 
Although it has been demonstrated that contributions 
from internal friction cause a sharp increase in C, for very 
short MTs, may be presumed to obey Eq.Q in the 
range of Lc considered in the current study [18 . We 
modelled damping according to the Rayleigh damping 
model: 



C(e) aM(e) 



•/3K 



(e), 



(8) 



where C(e) is the elements damping matrix, M(g) is the 
elements mass matrix and K(g) is the elements stiffness 
matrix. In our particular case of linear Timoshenko beam 
elements, we set a = C/p where p is the density of the 
MT and A is the cross-sectional area, and /3 = 0. 



2. Solving 

Each realization was solved by the Abaqus/Standard 
FE solver on a Dell Precision T7500 workstation (Xeon 



X5680, two CPUs @ 3.33 GHz). The Abaqus/Standard 
FE solver is based on an efficient implicit time integration 
scheme, which makes it suitable for low-speed dynamic 
events, such as thermal undulations of polymers. The 
handling of the realizations (automatically starting the 
solver and collecting the data afterwards) is performed 
by a Python script. 

For solving the model, we introduced two different 
timeframes: the simulation time frame, tgim, and the data 
collection time frame tcoii- At tgim = 0, the simulation 
is initiated and the initially straight contour of the fila- 
ment starts to show bending fluctuations. To ensure that 
data collection is started after the polymer has reached 
equilibrium with its environment, we began data collec- 
tion after the largest relaxation time of the polymer, Tc, 
elapsed, at which instant we set tcow (Figure [4|. This 
largest relaxation time of the polymer is determined by 
i3il5j 



-^qi~^ = ^^^^ 



(9) 



where the wave vector of the first bending mode in 
Eq.([9| depends on the imposed boundary conditions. For 
example, for hinged boundary conditions = l.bir/L ^ 
[34 and for free boundary conditions ql = tt/L ^ [5^ (Fig- 
ure [3|. 
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FIG. 3. Two sets of boundary conditions have been used for 
validation: a freely floating polymer (top) and a polymer with 
hinged ends (bottom). 



time: t. =0 t =0 

Sim coll 



30 T 



stepsize: Atp^=T/103 M^=xjyQ>^ l^\=xm ^\=xJW 

FIG. 4. A schematic depiction of the simulation timeline in- 
cluding the simulation time frame and collection time frame. 
The step size is adaptive: after relaxation is performed with 
a coarse step, the framework switches to a much smaller step 
size, which is then gradually increased. 

For various purposes, the dynamics need to be well de- 
scribed at very small timescales (t « Tc) as well as very 
large timescales {t » Tc) (e.g. for plotting the MSD of 
the end-to-end distance). In order to optimize the reso- 
lution for all timescales, variable time stepping is imple- 
mented: the time step that is used, depends on Tc and 
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the advancement of the reahzation. The procedure is as 
fohows: first, the relaxation time is calculated based on 
the input parameters, Eq.([9|. A coarse time stepping 
resolution is used to let the polymer equilibrate with its 
fluid environment (0 < t < Tc). Once the equilibrium 
has been obtained, the FE solver switches to the high- 
est resolution and accurately captures the dynamics at 
the smallest timescale. The time resolution is gradually 
decreased, every decade of timescale, up to the largest 
time step that still allows for proper convergence of the 
realization (Figure [4|. 

In the current FE framework, the computational cost 
of the simulation is, due to the variable time stepping, 
to a great extent independent of the relaxation time. 
Therefore, the computational effort is largely indepen- 
dent of the contour length of the polymer, the boundary 
conditions, the drag coefficient and the bending rigidity. 
However, for very large relaxation times or very low poly- 
mer rigidity the gain of variable time stepping is limited, 
due to a limit to the maximum step size that still allows 
for proper FE convergence. As an indication, a typi- 
cal realization (SOtc of simulated time) of an isotropic 
polymer lasted ± 45 minutes on the specified computer 
system. For high levels of anisotropy or very long relax- 
ation times, the duration of a realization could take up 
to 4 hours. 

The realizations within each ensemble are independent 
of each other (due to the ergodic nature of the studied 
phenomena). Therefore, simulations can be parallelized, 
reducing the overall simulation time. 



3. Post-processing 

The solver outputs data files that contain a global time 
stamp of each step AtpE and the corresponding positions 
of the first and last nodes of the polymer. 

In order to quantify the time evolution of the poly- 
mer fluctuations, we define the mean square displacement 
(MSD) of the end-to-end distance as 

(SRHtcou)} = {[Ritcou) - R{tcoU = 0)]'), (10) 

where R{t coii) is the end-to-end distance or projected 
length of the filament (Figure [5| . As can be seen from 
Eq.([lO|, the MSD of the end-to-end distance is zero if 
^coii = and gradually increases with time. The end-to- 
end distance is determined by subtracting the position 
vector of the last node and the position vector of the 
first node of the polymer. This procedure is performed 
for all time steps AtpE- The end-to-end distance at the 
of the collection R{tcoi\ = 0) is subtracted from the total 
end-to-end distance and the resulting values are squared 
to obtain the displacement of the end-to-end distance 
SR'^itcow)- We then average over all Nr realizations in 
the ensemble to obtain the MSD of the end-to-end dis- 
tance {SR'^ {t coll))' 




FIG. 5. For the MSD of the end-to-end distance, that was 
used as observable in this framework, we collect the coordi- 
nates of the first and last node. For tcoii = (SR^ (tcoii)) is 
set to zero and increases as time evolves. 



B. Validation results 

For validation purposes, two sets of boundary condi- 
tions have been considered: a freely floating polymer and 
a polymer with hinged ends (Figure [3|. In case of hinged- 
hinged boundary conditions, the end-nodes of the poly- 
mer are allowed to move along the x-axis, but are con- 
strained to move along the y-direction. 

According to the theory of the Brownian dynamics of 
polymers [ 341436] two regimes can be identified in the 
time evolution of the mean square end-to-end distance: 
for very small times t « Tc [37l, (5i?^(t)) increases, 
obeying a power law that scales with t^^^ , whereas 
for t » (^5R^{t)) reaches a universal equilibrium 
value of {SR^{t))^^ = 1^^/901^ We rescaled the 

simulated (^6R'^{t)) and analytical solution by defining 
F{t) = (SR'^it)) [ml/Lc^) and we rescaled time by 
defining t = t/rc. For both boundary conditions (free 
and hinged) and contour lengths (10 pm and 20 pm) a 
good agreement with the theoretical predictions of the 
WLC model [3T, ""36 is observed for over 5 decades of 
timescale. The jagged shape of the simulated curves is 
due to the random nature of the observed phenomenon 
and is expected to die out completely for a higher number 
of realizations, as shown in a similar research [23]. It is 
not the goal of this particular study to show convergence 
of the finite element towards the analytical solution for 
a large number of realizations (this has been done before 
in [23 for = 4000). For the purpose of the cur- 
rent investigation, namely validation of the framework, 
Nji = ±50 realizations per ensemble was considered suf- 
ficient to show a good correspondence with the WLC 
theory. 



C. Anisotropic polymers 

For the validation of the framework we presumed 
isotropic material properties. However, due to their 
atomic architecture (Figure [!}, MTs are known to be 
anisotropic: the longitudinal Young's modulus is deter- 
mined by the strength of the head-tail a/S — a/S tubulin 
bonds, whereas the shear modulus is determined by the 
much weaker inter-protofilament bonds |14Hl6j . Further- 
more, in the plane of the transverse cross- sect ion, MTs 
give an isotropic response to stress. From here on, we 
introduce a transversely isotropic material model, which 
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FIG. 6. Rescaled MSD of the end-to-end distance, F{t) = 
{SR^{tcon)) (90-£p/Lc'^) , as a function of rescaled time t/rc- 
The plots show the F{t) for various combinations of boundary 
conditions and contour lengths [(a) BC: free, Lc = 10 pm, (b) 
BC: hinged, Lc = 20iim]. 



accounts for the anisotropic material behaviour along the 
backbone of the MTs and isotropic properties along their 
cross-section. This approach is similar to other studies 
on MTs |2l [28] . According to linear elasticity theory, the 
elastic compliance matrix of a transversely isotropic ma- 
terial is defined by the Young's moduli in the plane of 
isotropy, E2 = Es = Ep^ the transverse Young's modulus 
El = Ef ^ the Poisson ratio's u 



'p 1 ^pt 



, Utp and the in- 
plane and transverse shear moduli, Gis = G23 = Gp and 
G12 = Gt, respectively. Therefore, the 3-D stress-strain 
law reduces to: 



1/Et -i^tp/Et -i^tp/Et 
-Vpt/Ep l/Ep —Vp/Ep 
-Vpi/Ep —Up/ El 1/Ep 













[cr 33) 







713 




\723y 





1/Gt 








(12) 



Because of the 2-D confinement of the simulation, deflec- 
tion of the filament is only determined by the longitudi- 
nal elastic modulus, Ei, and longitudinal shear modulus, 
G12. Therefore, we introduced the transverse isotropy 
through the parameter, which sets the order of the 
ratio between the longitudinal shear modulus and corre- 
sponding Young's modulus. 



G12 



10^. 



(13) 



To conform our research as much as possible to the exper- 
imental work that has been done on grafted microtubules 
[2] , we clamped the left-end of the MT by applying con- 
straints in all three degrees of freedom (2 lateral, 1 rota- 
tional), whereas we left the right-end unconstrained (Fig- 
ure [7|). This new set of boundary conditions changes the 
longest wave vector, which is now given by ql = 1.875/Lc 
j39] . Substitution of this wave vector in Eq.([9| gives us 
for an isotropic grafted polymer a relaxation time of 



EI 



1.875 ) 



(14) 



For anisotropic polymers with a low shear modulus, it is 
expected that the effective bending rigidity decreases and 
thus the relaxation time increases. To estimate the effec- 
tive relaxation time of an anisotropic polymer, r^^, we 
studied the MSD of the end-to-end distance, (^i?^(tcoii)) 5 
for each unique set of parameters. Based on the onset 
of the {^R^{tco\\)) plateau we estimated rf^. For sub- 
sequent simulations, data was only collected after this 
effective relaxation time plus an additional safety mar- 
gin elapsed. We then considered for each ensemble the 
position of the grafted MT tip in the x-y plane which we 
sampled at an interval of This allowed for ac- 

curate measurement of the distribution function P(x, y) 
and transverse displacement ^^of the tip position (Fig- 
ure [t]). Then, by integrating the probability distribution 
function over the y-axis, we found the asymmetric prob- 
ability density function along the x-axis P{x). 




(11) 



FIG. 7. Schematic drawing of the studied microtubule: the 
left-end is constrained in all degrees of freedom, whereas the 
right-end remains unconstrained. The mean square of the 
transverse displacement of the tip, is used as a measure 
for the persistence length. 

For semiflexible polymers in the stiff limit (Lc < 
^p), this asymmetric distribution is peaked toward full 
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stretching and has a typical width of [iU] : 



P(x,y) 



(15) 



Similarly, by integration along the x-axis, we obtained 
the distribution function P(^), which is a Gaussian dis- 
tribution centred at = 0, the variance of which is given 
by the mean square transverse displacement {y\) [40] . 
This mean square transverse displacement {y\) can be 
directly related to the persistence length of the polymer, 
according to [5l [4Ql[4T] : 



Li 



(16) 



For each realization that is added to the ensemble, the 
mean square transverse displacement, and thus the per- 
sistence length, changes. The ensemble was considered 
to contain enough realizations for a reliable estimate of 
the persistence length, if for each consecutive addition of 
the last five realization to the ensemble the mean square 
transverse displacement did not change by more than 3%. 



III. RESULTS 

First, to affirm our approach, we sampled the tip 
of grafted isotropic MTs of various lengths and calcu- 
lated the simulated persistence length using Eq.(16). We 
found a good agreement with the theoretical persistence 
length (^p = EI/ksT = 6.3mm). As an illustration. 
Figure [8] shows the tip position of an isotropic MT of 
contour length Lc = 10 pm in the x-y plane, based on 
1.8 X 10^ samples. By integration along the y-axis, we 
find the characteristic non-Gaussian and asymmetric dis- 
tribution function Pfx) with a width of Ly = 0.015 pm, 
as predicted by Eq.([l5| according to the WLC theory 
[31]. Similarly, integration along the x-axis results in 
the Gaussian distribution P(y), for which a good agree- 
ment is observed with the WLC theory (Figure [8|. In- 
deed, we find a mean square transverse displacement of 
(yl) = 0.0550 pm, which is within 4% accuracy of the 
theoretical value of (?/^) = 0.0529. 

We calculated the effective persistence length for MTs 
of various contour lengths and orders of anisotropy by 
means of Eq. ( 16 ) . These results are presented in Figure|9] 
We find that with increasing orders of transverse isotropy 
(x = 1,2,3) the persistence length decreases for all con- 
tour lengths. In this range of anisotropy, short MTs seem 
to be less affected by a lower longitudinal shear modulus 
than long MTs. For highly anisotropic MTs (x = 6), 
we see a drop in persistence length of two orders com- 
pared to their isotropic counterpart. Although changes 
in persistence length with contour length are seen for 
all values of x? this dependency is monotonic only for 
X = 6. For this value, a gentle, but noticeable, increase 
of the persistence length with increasing contour length 



is observed. In Figure 10, the contour length-dependency 
of the persistence length of MTs with the highest level 
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FIG. 8. a) Point cloud of sampled tip position, based on 
1.8 X 10^ samples, b and c) After integration along the y and 
X-axis we find the two characteristic distributions P{x) and 
P{y), respectively, c) The solid line shows the analytical pre- 
diction according to the worm-like chain model. 



of anisotropy (x = 6) of the FE simulations is compared 
to the simplified cantilever beam model, based on the 
Timoshenko beam formalism and single static tip load, 
as used in other studies [2j [18] . 



IV. DISCUSSION 

Certain biopolymers, such as MTs, display an 
anisotropic material response, which is inherent to their 
molecular structure O [5]. Unfortunately, the effective 
persistence length of such anisotropic filaments cannot 
be explicitly expressed in terms of their material prop- 
erties. Additionally, the anisotropic material properties 
and relative thickness of MTs cannot be accurately im- 
plemented in conventional numerical and analytical mod- 
els, since such models are inherently built on slender- 
ness and isotropy approximations. Therefore, researchers 
have been compelled to simple static models to interpret 
their experimental observations. For example, based on 
such a static model, unexplained length-dependency of 
the persistence length has been related to anomalously 
high levels of anisotropy [2 . However, the validity of 
such static models, which are based on small angle ap- 
proximations and a single deterministic tip force, may 
be questionable when describing the finite displacements 
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FIG. 10. Persistence length versus contour length for MTs 
with a six-order lower longitudinal shear modulus than cor- 
responding Young's modulus according to the static Tim- 
oshenko beam model ^ |18] with single tip force (dashed) 
and dynamic FE simulations with distributed random forces 
(solid). 



that are typically observed for highly anisotropic fila- 
ments under distributed random forces. To overcome 
these limitations we developed and validated a frame- 
work, based on the FE method set out by Cyron and Wall 
in ^23^. This FE framework accounts for the interaction 



between semifiexible polymers and solvent molecules by 
random force generation and application. The FE model 
also allows for the implementation of advanced mate- 
rial models and nonlinear elasticity. Additionally, the 
developed framework is computationally efficient com- 
pared to conventionally used bead models. By apply- 
ing the FE framework to isotropic MTs of various con- 
tour lengths, we found that the time evolution of thermal 
fluctuations is in good agreement with the predictions of 
the WLC model [H |36]. Furthermore, the FE frame- 
work allowed us to implement various levels of mate- 
rial anisotropy. By studying the distribution function of 
the tip of grafted isotropic and anisotropic MTs, we cal- 
culated the persistence length for four different contour 
lengths, {Lc = 5, 10, 15,20 pm). For isotropic MTs with 
a contour length of Lc = 10 pm, we found a good corre- 
spondence with theoretical predictions and other Monte 
Carlo studies |4Q] . 

Three trends were identified from the simulation re- 
sults. First, the implementation of a transversely 
isotropic material model with a lower longitudinal shear 
modulus than corresponding Young's modulus, results 
in a decrease of the persistence length for all contour 
lengths. Second, for the first three orders of transverse 
isotropy (x = 1, 2, 3), this reduction in persistence length 
is greater for long MTs than for short MTs. Third, imple- 
menting the highest order of transverse isotropy (x = 6) 
yields a two order decrease in persistence length, as com- 
pared to the isotropic case. Additionally, for x = 6, a 
slight increase of the persistence length with increasing 
contour length is observed. 

The two orders lower persistence length of highly 
anisotropic MTs (x = 6) is of the same order as the 
persistence length reported in experimental studies for 
short and very short MTs |^. These experimen- 
tal results have been explained by considering MTs as 
weakly- coupled (at intermediate length scales [18 ) or de- 
coupled (at short length scales [42 ) assemblies of protofil- 
aments. We can thus confirm that accounting for inter- 
protofilament decoupling by lowering the longitudinal 
shear modulus indeed results in such short persistence 
lengths. However, the above-mentioned experiments also 
confirmed a significantly longer persistence length of 4 — 8 
mm for long MTs of Lc = 20 — 40 mm, as also estab- 
lished in earlier studies, e.g. [3 . Such a striking in- 
crease in persistence length has been attributed to the 
presumption that, in this length regime, MTs behave 
as fully coupled homogenous structures with negligible 
inter-protofilament sliding. From our simulation results, 
we conclude that a six order difference in longitudinal 
shear modulus and corresponding Young's modulus alone 
cannot account for this reported length-dependence of 
the persistence length. Our findings agree well with pre- 
vious research on MT rigidity, based on experimentally 
measurement of buckling force [33, 34] and mode decom- 
position of free-floating MTs [3 , in which no significant 
length-dependency of flexural rigidity was reported. Fur- 
thermore, the outcomes of the current study reveal the 
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limitations of the simplified mechanical cantilever model 
that was used in other studies |2l [18] to interpret the 
observed length- dependency. Such a static Timoshenko 
beam model is based on small angel approximations and 
a single tip load. Therefore, it is unsuitable for accu- 
rately capturing the intricate dynamics and large dis- 
placements of a highly anisotropic MTs in equilibrium 
with their fluid environment. The current research sheds 
new light on the unexplained discrepancy between the 
indirectly deduced longitudinal shear modulus (six order 
difference with corresponding Young's modulus), based 
on such a simplified mechanical model, and direct exper- 
imental and computational measurements of the longitu- 
dinal shear modulus [l3l[T7] (up to three orders difference 
with corresponding Young's modulus). The outcomes of 
the current study should be interpreted as an encour- 
agement to consider other causes for the experimentally 
observed length-dependence than high anisotropy alone. 

The current study is subject to several limitations. 
A first limitation is that internal friction due to liquid 
flowing through narrow pores of the MT, as pointed 
out in [18], was neglected. Including internal friction 
would have caused a sharp peak in the drag coefficient 
for very short MTs. However, in the range of contour 
lengths considered in the current study, the drag coeffi- 
cient may be presumed constant [18 . A second limitation 
of the current study is that we approximated friction by 
an isotropic friction model. Although an isotropic fric- 
tion model serves as a good first approximation, [23 an 
anisotropic friction model, that takes into account differ- 
ent friction coefficients perpendicular and longitudinal to 



the filament, enhances the accuracy of the absolute val- 
ues of the model, see [24- Additionally, we only implic- 
itly accounted for the hydrodynamic interactions between 
distant polymer sections through scaling of the drag coef- 
ficient by a logarithmic length-dependent correction fac- 
tor. Indeed, we noticed small differences in the amplitude 
of the saturation plateau of the MSD of the end-to-end 
distance upon changing the drag coefficient correction 
factor. Such an approximation is common practice for 
stiff polymers, because for such stiff polymers remote in- 
teractions, such as one segment shielding another seg- 
ment, will be negligible [431I11]- However, for long and 
flexible polymers, hydrodynamic interaction between seg- 
ments as well as self-avoidance may become appreciable 
and should be taken into account. In the current re- 
search the most flexible polymer {Lc = 20 pm, X = 6), 
was still in the stiff regime, ^ < 0.1. Therefore, it is 
expected that neglecting hydrodynamic interaction and 
self- avoidance will only have had a minor effect on the 
simulated values and a negligible effect on the observed 
trends [40 . Finally, due to high computational costs, the 
number of realizations that were performed per ensemble 
was limited. A higher number of realizations would have 
further improved the accuracy. 

Despite the above mentioned limitations, the accuracy 
of the model exceeds largely the required accuracy to 
support our main finding: a high anisotropy cannot ac- 
count for the large length-dependence of the persistence 
length. 
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